Genetic algorithms and supernovae type la analysis 
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We introduce genetic algorithms as a means to analyze supernovae type la data and extract 
model-independent constraints on the evolution of the Dark Energy equation of state w(z) = . 
Specifically, we will give a brief introduction to the genetic algorithms along with some simple 
examples to illustrate their advantages and finally we will apply them to the supernovae type la 
data. We find that genetic algorithms can lead to results in line with already established parametric 
and non-parametric reconstruction methods and could be used as a complementary way of treating 
SNIa data. As a non-parametric method, genetic algorithms provide a model-independent way to 
analyze data and can minimize bias due to premature choice of a dark energy model. 
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I. INTRODUCTION 

The accelerated expansion of the universe has been 
confirmed during the last decade by several observational 
probes [11, IS Mi, [3 . The origin of this acceleration 
may be attributed to either dark energy with negative 
pressure, or to a modification of General Relativity that 
makes gravity repulsive at recent times on cosmological 
scales. One way to distinguish between these two possi- 
bilities and identify in detail the gravitational properties 
of dark energy or modified gravity is to have a detailed 
mapping of the expansion rate H{z) as a function of the 
redshift z. This is equivalent to identifying the dark en- 
ergy equation of state w{z) = ^ which can be written 
as 

w{z) = -1 + -(1 + z) , (1.1) 

6 amz 

where SH{z)^ = H (z)"^ / Hq - flomi'^ + z)^ accounts for all 
terms in the Friedmann equation not related to matter. 
The cosmological constant {w{z) = —1) corresponds to a 
constant dark energy density, while in general it can be 
time dependent. 

It has been shown [1] that a, w{z) observed to cross 
the line w{z) = — 1 (phantom divide line) is very hard 
to accommodate in a consistent theory in the context of 
General Relativity. On the other hand, such a crossing 
can be easily accommodated in the context of extensions 
of General Relativity 0. Therefore, the crossing of the 
phantom divide line w = —1 could be interpreted as a 
hint in the direction of modified gravity, see for example 
0, Such a hint would clearly need to be verified 

by observations of linear density perturbation evolution 
through e.g. weak Icnsing 12] or the redshift distortion 
factor 131. 
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Early SNIa data put together with more recent such 
data through the Gold dataset [J, [l^ have been used 
to reconstruct w{z), and have demonstrated a mild pref- 
erence for a w{z) that crossed the phantom divide line 
(Til [TH - A cosmological constant remained consistent 
but only at the 2a level. However, the Gold dataset has 
been shown to suffer from systematics due to the inho- 
mogeneous origin of the data [lB|. More recent SNIa 
data (SNLSpJ], ESSENCE[i|, HSTfli]) re-compiled in 
[l9l | have demonstrated a higher level of consistency with 
ACDM and showed no trend for a redshift dependent 
equation of state. 

However, despite the recent progress the true nature 
of Dark Energy still remains a mystery with many pos- 
sible candidates being investigated, see for example 20l |. 
The simplest possible candidate assumes the existence of 
a positive cosmological constant which is small enough 
to have started dominating the universe at recent times. 
This model provides an excellent fit to the cosmological 
observational data [2l[ and has the additional bonus of 
simplicity and a single free parameter. Despite its sim- 
plicity and good fit to the data, this model fails to explain 
why the cosmological constant is so unnaturally small as 
to come to dominate the universe at recent cosmological 
times, a problem known as the coincidence problem and 
there are specific cosmological observations which differ 
from its predictions (2^.(23|. 

A further complication in the investigation of the be- 
havior of dark energy occurs due to the bias introduced 
by the parameterizations used. At the moment, there is a 
multitude of available phenomenological ansatze for the 
dark energy equation of state parameter w, each with 
its own merits and limitations (see [l^l and references 
therein). The interpretation of the SNIa data has been 
shown to depend greatly on the type of parametrization 
used to perform a data fit [ISl ■ Choosing a priori a model 
for dark energy can thus adversely affect the validity of 
the fitting method and lead to compromised or mislead- 
ing results. The need to counteract this problem paved 
the way for the consideration of a complementary set 
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of non-parametric reconstruction techniques [1^, [2^, [l^l ■ 
These try to minimize the ambiguity due to a possibly 
biased assumption for w by fitting the original datasets 
without using any parameters related to some specific 
model. The result of these methods can then be inter- 
preted in the context of a DE model of choice. Non- 
parametric reconstructions can thus corroborate para- 
metric methods and provide more credibility. However, 
they too suffer from a different set of problems, mainly 
the need to resort to differentiation of noisy data, which 
can itself introduce great errors. 

In this paper, we consider a method for non-parametric 
reconstruction of the dark energy equation of state pa- 
rameter w, based on the notions of genetic algorithms 
and grammatical evolution. We show how the machin- 
ery of genetic algorithms can be used to analyze SNIa 
data and provide examples of its application to deter- 
mine dark energy parameters for specific models. Our 
results are consistent with already known findings of both 
parametric and non-parametric treatments. In particu- 
lar, we obtain parameters which agree within la with 
AC DM, but show a slight trend towards a varying pa- 
rameter w{z). As a non-parametric method, genetic al- 
gorithms provide a model-independent way to analyze 
data and can minimize bias due to premature choice of a 
dark energy model. 

A rather formal definition of the Genetic Algorithms 
(GA) [2^ is that of a computer simulation in which a 
solution to some predefined problem is searched for in 
a stochastic way by evolving an ensemble (population) 
of candidate solutions (usually called chromosomes) to- 
wards an optimal state, where some of the chromosomes 
have reached sufficiently close to the true solution. This 
approach to problem solving resembles the evolutionary 
pattern encountered in nature, where an optimal or dom- 
inant organism emerges through a series of apparently 
random mutations and combinations between different 
individuals. The concept of optimality of a chromosome 
is usually encoded in the genetic algorithm language in a 
fitness function, corresponding to the phenotype of an 
individual. Candidate solutions with a higher fitness 
function are considered to be better than the rest of the 
population, in the sense that they are closer to the true 
solution. Genetic operators such as mutation and combi- 
nation (crossover) are then implemented to produce off- 
spring from a given population with better properties. 
Chromosomes with a good fitness are usually more likely 
to crossover and produce descendants and the combina- 
tion often takes place between the best candidate solu- 
tions of a generation. It is thus more likely for descen- 
dants to bear a combination of already favorable genes, 
leading to an improvement in fitness. The genetic algo- 
rithm method is in a way similar to search algorithms, 
such as hill climbing, A* [2^ or simulated annealing's^. 

GAs are more useful and efficient than usual techniques 
when 

• The parameter space is very large, too complex or 
not enough understood. 



• Domain knowledge is scarce or expert knowledge is 
difficult to encode to narrow the search space. 

• Traditional search methods give poor results or 
completely fail. 

so naturally they have been used with success in many 
fields where one of the above situations is encountered, 
like the computational science, engineering and eco- 
nomics. Recently, they have also been applied to study 
high energy physics [3l|,[33|, HI] gravitational wave de- 
tection 34] and gravitational lensing [3^. Since the na- 
ture of Dark Energy still remains a mystery, this makes 
it for us an ideal candidate to use the GAs as a means 
to analyze the SNIa data and extract model independent 
constraints on the evolution of the Dark Energy equation 
of state w{z). In Section 2 we will briefly describe the 
standard way of analyzing SNIa data while a more thor- 
ough treatment of GAs and grammatical evolution (GE) 
is deferred to Section 3. Section 4 contains the general 
methodology of the GA paradigm and its application to 
Supernovae data. 

II. THE STANDARD SNIa DATA ANALYSIS 

We use the SNIa dataset of Kowalski et. al. [4l| con- 
sisting of 414 SNIa, reduced to 307 points after various 
selection cuts were applied in order to create a homoge- 
neous and high-signal-to-noise dataset. These observa- 
tions provide the apparent magnitude m(z) of the super- 
novae at peak brightness after implementing the correc- 
tion for galactic extinction, the K-correction and the light 
curve width-luminosity correction. The resulting appar- 
ent magnitude m{z) is related to the luminosity distance 
Dl{z) through 

mtUz) = M{M, Ho) + 5logw{DLiz)) , (2.1) 

where in a flat cosmological model 

DLiz) = (1 + z) /' dz'——^ ^ , (2.2) 

Jo H(z';iloni,wo,wi) 

is the Hubble free luminosity distance {Hod if), and M is 
the magnitude zero point offset and depends on the abso- 
lute magnitude M and on the present Hubble parameter 
Ho as 

M ^ M + 5/0510 (^)+ 25 = 

= A/ - 5/0510/1 + 42.38. (2.3) 

The parameter M is the absolute magnitude which is 
assumed to be constant after the above mentioned cor- 
rections have been implemented in m{z). 

The SNIa datapoints are given, after the corrections 
have been implemented, in terms of the distance modulus 

Hobs{zi) = mobs{zi) ~ M. (2.4) 
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The theoretical model parameters are determined by 
minimizing the quantity 

XSNIa{^Om,Wo,Wi) = 2 , (2.5) 

where N = 307 and ct^ ^ are the errors due to flux uncer- 
tainties, intrinsic dispersion of SNIa absolute magnitude 
and peculiar velocity dispersion. These errors are as- 
sumed to be Gaussian and uncorrelated. The theoretical 
distance modulus is defined as 

lith{zi) = mth{zi) - M = 5logio{DL{z)) + , (2.6) 

where 

/xo = 42.38 - 5Zogio/i , (2.7) 

and fiobs is given by (|2.4p . The steps we followed for the 
usual minimization of (|2.5p in terms of its parameters are 
described in detail in Rcfs. [ll[ll|4l|. 

This approach assumes that there is some theoreti- 
cal model available, given in the form of H{z]pi) and in 
terms of some parameters pi like the f2om or the wq and 
wi of the CPL ansatz, which is to be compared against 
the data. As a result of the analysis, the best-fit pa- 
rameter values and the corresponding la error bars are 
obtained. The problem with this method is that the ob- 
tained values of the parameters are model-dependent and 
in general models with more parameters tend to give bet- 
ter fits to the data. This is where the GA approach starts 
to depart from the ordinary parametric method. Our 
goal is to minimize a function like ()2.5|) . not using a can- 
didate model function for the distance modulus and vary- 
ing parameters, but through a stochastic process based 
on a GA evolution. This way, no prior knowledge of a 
dark energy models is needed to obtain a solution. The 
resulting "theoretical" expression for fithiz) is completely 
parameter-free. This is the main reason for the use of GA 
in this paper. We now proceed to give a more thorough 
treatment of GAs and the results of their application in 
the next sections. 



III. GENETIC ALGORITHMS 

A. Overview 

GAs were introduced as a computational analogy of 
adaptive systems. They are modelled loosely on the prin- 
ciples of the evolution via natural selection, employing a 
population of individuals that undergo selection in the 
presence of variation-inducing operators such as muta- 
tion and crossover. The encoding of the chromosomes is 
called the genome or genotype and is composed, in loose 
correspondence to actual DNA genomes, by a series of 
representative "genes" . Depending on the problem, the 
genome can be a series (vector) of binary numbers, deci- 
mal integers, machine precision integers or reals or more 



complex data structures. As mentioned in the introduc- 
tion, a fitness function is used to evaluate individual chro- 
mosomes and reproductive success usually varies with fit- 
ness. The fitness function is a map between the gene se- 
quence of the chromosomes (genotype) and a number of 
attributes (phenotype) , directly related to the properties 
of a wanted solution. Often the fitness function is used to 
determine a "distance" of a candidate solution from the 
true one. Various distance measures can be used for this 
purpose (Euclidean distance, Manhattan, Mahalanobis 
etc.). 

The algorithm begins with an initial population of can- 
didate solutions, which is usually randomly generated. 
Although GA's are relatively insensitive to initial condi- 
tions, i.e. the population we start from is not very sig- 
nificant, using some prescription for producing this seed 
generation can affect the speed of convergence. In each 
successive step, the fitness functions for the chromosomes 
of the population are evaluated and a number of genetic 
operators (mutation and crossover) are applied to pro- 
duce the next generation. This process continues until 
some termination criteria is reached, e.g. obtain a solu- 
tion with fitness greater than some predefined threshold 
or reach a maximum number of generations. The later is 
imposed as a condition to ensure that the algorithm ter- 
minates even if we cannot get the desired level of fitness. 

The various steps of the algorithm can be summarized 
as follows: 

1. Randomly generate an initial population Af(0) 

2. Compute and save the fitness for each individual m 
in the current population AI{t). 

3. Define selection probabilities p{m) for each individ- 
ual m in M{t) so that p{m) is proportional to the 
fitness. 

4. Generate M{t+ 1) by probabilistically selecting in- 
dividuals from M(t) to produce offspring via ge- 
netic operators. 

5. Repeat step 2 until satisfying solution is obtained, 
or maximum number of generations reached. 

The paradigm of GAs descibed above is usually the one 
applied to solving most of the problems presented to GAs. 
Though it might not find the best solution, more often 
than not, it would come up with a partially optimal so- 
lution. 



B. Grammatical evolution 

The method of grammatical evolution (GE) was 
first introduced as an alternative approach to genetic 
programming [3^ [s^- Genetic programming is a varia- 
tion of GAs, which uses programming trees as chromo- 
somes instead of vectors of binary numbers or strings. 
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These chromosomes can be used to express formal expres- 
sions, mathematical functions or even entire programs 
written in some specific programming language. GA op- 
erators such as mutation or crossover can then be inter- 
preted as operations which alter or exchange subtrees of 
a chromosome tree, or splice together two distinct trees. 

Grammatical evolution can also be used to produce ar- 
bitrary programs in any programming language but with- 
out the need for tree structures, which usually require 
complicated and cumbersome methods to implement ge- 
netic operators. Instead, one can use GE to produce an 
expression or program using an ordinary binary vector 
chromosome and a generative grammar, which maps the 
chromosome to an expression using a set of grammatic 
rules. These rules are usually expressed in a Backus-Naur 
form (BNF), which can be used to recursively convert a 
binary vector into a string of terminal symbols, which 
corresponds to the output. In this way, GE introduces 
an additional layer of abstraction between the genome 
(vector of binaries) and the phenotype (fitness). Usu- 
ally, the transition between the two is accomplished in 
one step, by evaluating the fitness function of a chro- 
mosome. In the case of grammatical evolution, there is 
the additional step of converting the initial chromosome 
into an expression using the available grammar. The ex- 
pression is then evaluated to find the fitness of the chro- 
mosome. GA operators can be easily implemented in 
this scheme on the level of vector chromosomes, since 
the formal expressions are evaluated in the next step and 
are not necessary to perform the mutation and crossover. 
The added abstraction layer also is more similar to what 
occurs in nature, where a transcription process is used 
to map DNA base sequences (genotype) into proteins 
(phenotype) . GE has already been applied to problem 
such as symbolic regression [37|, minimization[40], robot 
control [3811 and finance 1391. 



C. Initialization 

The first stage of the algorithm involves the creation 
of an initial population. The chromosomes are generated 
in a random way and in large enough quantities. The 
size of the population is very important to successfully 
mimic a natural selection process, as diversity is crucial 
to produce a variety of genes and phenotypes. It is thus 
advisable to use as large a population sample as possible 
in a simulation, while keeping computational complexity 
under control. Usually the chromosome count remains 
fixed throughout the evolutionary phase. Typically this 
involves a few hundreds, up to several thousands of indi- 
viduals, although these numbers vary with the type of the 
problem and the resources available. Similarly, the size 
of each chromosome (number of genes) is going to affect 
the speed of the algorithm and the quality of the final 
solution. Large genomes lead to richer phenotypes and 
better chances to increase the fitness through mutation 
and crossover, at the expense of computational efficiency. 



D. Selection 

The purpose of this phase is to select the chromosomes 
in the population which are going to breed to give a new 
generation. This is the step where the equivalent of nat- 
ural selection takes place. Chromosomes are chosen for 
crossover based on their fitness. Solutions with a better 
fitness function are the dominant members of the popula- 
tion, hence they are more likely to be selected to produce 
offspring. 

Selection is performed in a stochastic way to avoid pre- 
mature convergence on suboptimal solutions. The algo- 
rithm assigns some selection probability to each chromo- 
some according to its fitness function. The totality of 
the population or some part of it can be used to draw 
a number of candidates for crossover. The most popular 
and heavily used methods for selection are: 

• Roulette wheel selection. 

• Tournament selection. 

In the case of roulette wheel selection, a probability for 
crossover of a chromosome is considered to be propor- 
tional to its fitness. This is a widely used method, but 
suffers from a number of drawbacks. Since chromosomes 
with higher fitness are more likely to be selected for 
crossover, once a suboptimal solution occurs, it may dom- 
inate the mating pool. At the same time, less favored 
chromosomes tend to be completely neglected as candi- 
dates for breeding. This leads to premature convergence 
to suboptimal solutions and stagnation of the evolution- 
ary computation. Roulette wheel selection is also not so 
easy to implement. 

On the other hand, tournament selection chooses the 
candidates for crossover by taking smaller groups of ran- 
domly selected individuals and picking the dominant 
member of each group. This is an easy algorithm to 
implement and also more in line with what we usually 
encounter in nature. Suboptimal chromosomes are thus 
more likely to be picked than in the roulette wheel case 
and premature convergence can be avoided more effec- 
tively. Tournament selection is going to be used in our 
treatment of SNIa data. 

E. Reproduction 

After selecting the candidates for breeding, the next 
step is to combine the corresponding genomes to produce 
the members of the next generation. This is achieved 
by the application of genetic operators, mutation and 
crossover. 

For each pair of parent chromosomes picked from the 
population, two children chromosomes are produced by 
exchanging the genomes of the two parents. The eas- 
iest and more popular method is one-point crossover, 
where a particular point is chosen on the genome se- 
quence and the genes from that point on completely ex- 
changed between the chromosomes. The descendants 
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are then subjected to mutation with some predefined 
probabihty (mutation rate). This process of selec- 
tion/crossover/mutation is then repeated until enough 
individuals are produced to form a new generation. More 
often the entire population of the previous generation is 
replaced by members produced during the reproduction 
stage. In some cases, elitism is applied to ensure that 
some of the already found solution with the best fitnesses 
are passed on to the next generation. Although the two 
parent scheme is more biology inspired, recent researches 
suggested more than two parents may lead to better qual- 
ity chromosomes. 

By crossover and mutation, a new generation with dif- 
ferent members from the progenitor is created. New chro- 
mosomes tend to inherent the best qualities of the pre- 
vious generation, since members with higher fitness are 
more likely to be selected for reproduction. In general, 
the new chromosomes will have better fitnesses on the 
average and a better maximum fitness within the popu- 
lation. 



F. Termination 

The GA evolution continues until a termination con- 
dition has been reached. Some of the most common con- 
ditions for termination include: 

• An optimal solution is found based on some best 
fitness threshold. 

• The maximum number of generations is reached. 

• Allocated budget (CPU time/money) reached. 

• Further execution doesn't significantly improve the 
best fitness in the population. 

• Combinations of the above. 



IV. RESULTS 
A. General Methodology 

We first outline the course of action we follow to ap- 
ply the GA paradigm in the case of SNIa data. For the 
application of GA and GE on the dataset, we use a mod- 
ified version of the GDF^ tool ^, which uses GE as a 
method to fit datasets of arbitrary size and dimensional- 
ity. This program uses the tournament selection method 
for crossover. CDF requires a set of input data (train 
set), an (optional) test sample and a grammar to be used 
for the generation of functional expressions. The output 



^ |http://cpc.cs. qub .ac.uk/ summaries / ADXC] 



is an expression for the function which best fits the train 
set data. 

In our case, the train set is the SNIa dataset, con- 
taining 307 pairs of the distance modulus fi{z) and the 
redshift z as described in Section 2. We also modify GDF 
to read an additional file containing the errors of the cor- 
responding data points, which are used to evaluate the 
fitness function. No test dataset is needed for our pur- 
poses. 

The grammar used may involve standard numerical 
operations (-1-,—,*,/, etc.), as well as typical functions 
such as sin, cos, log, exp etc. and more complex produc- 
tion rules to turn non-terminal expressions into terminal 
symbols. The choice of the grammar vocabulary is ap- 
parently crucial for the behavior of the algorithm, the 
speed of convergence and the quality of the final fit. Our 
goal is to obtain a smooth fit for the distance modulus 
curve in the range of < z < 1.75, which will not lead to 
problematic behavior upon differentiation and can pro- 
duce a high goodness of fit. The use of the abs function 
is initially excluded, since as it turns out is prone to lead 
to discontinuities and unphysical behavior of the result- 
ing dark energy profile. Periodic functions may also lead 
to high frequency oscillatory behavior in the candidate 
solutions, as the fitting function attempts to interpolate 
between densely spaced fluctuating datapoints. Solutions 
of this sort generally seem to yield a better goodness of 
fit compared to grammars not employing periodic func- 
tions. Such oscillatory behavior is hard to distinguish 
within the accuracy available and there is little theoreti- 
cal justification to take it into account. We also find that 
this family of grammars in the end yields higher errors 
for the final model parameters, once a particular dark 
energy model is tested using the periodic solutions. We 
thus choose to exclude functions such as sin and cos for 
our final treatment of SNIa data. 

It is important to note that the method we discuss is 
highly non-parametric. By this we mean that no param- 
eters are used in the process of fitting the initial dataset 
to obtain a candidate solution. The GA does not apply 
some kind of minimization using a seed parametric func- 
tion. Instead, the maximization of the best population 
fitness, which is equivalent to a minimization proce- 
dure, is accomplished in a purely stochastic way. The 
resulting solution is given as an ordinary formal mathe- 
matical expression of a function of z, in which no param- 
eters appear. This should be contrasted to other non- 
parametric dark energy reconstruction methods, where 
some sort of parameter minimization takes place from 
the very beginning, even if these parameters are not di- 
rectly related to some specific dark energy model and 
only later do they get mapped to some corresponding set 
of model-dependent quantities [25]. This is an advantage 
of the GA approach, as it minimizes the amount of ini- 
tial model-dependent bias. Strictly speaking, the only 
parameter which can be tuned beforehand is the vocabu- 
lary of the grammar used. Grammatical evolution ame- 
liorates this potential source of bias, since it possesses 
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FIG. 1: Convergence of the genetic algorithm as a function of 
the number of generations. Curves are presented for different 
values of the total number of chromosomes in the population. 
The case with 190 chromosomes converges very slowly for the 
original SNIa dataset and is omitted. 

a far greater expressive power than ordinary fit metli- 
ods used, such as bin-wise polynomial fits. The apparent 
drawback of the GA is that it is in the end overly non- 
parametric, to the point where no established tools used 
in usual minimization techniques can be employed to es- 
timate and propagate errors. We will discuss later how 
this problem can be bypassed. The fitness function for 
the GA is chosen to be equal to — x|iv/a' where 

XSNIa — —2 ' l^-^J 

i=l % i 

and jl{z) — ^{z) — ^q. Notice the differences compared 
with (|2.5p . The XsNia doesn't depend on parameters 
and P,ga{z) is the reduced distance modulus obtained for 
each chromosome of the population through GE. We opt 
to fit for p.{z) instead of n{z), since the non-parametric 
approach we use prevents us from marginalizing over /iq. 

Another advantage is the fact that the GA does not 
require an assumption of fiatness as it can be seen from 
eq. (|2.5|) . since it was used in order to find the re- 
duced distance modulus Flatness (or non-flatness) 
comes into play when one tries to find the luminosity dis- 
tance D]^ (z) and the underlying dark energy model given 
fJ.GA{z). So, in this aspect the GA is independent of the 
assumption of flatness as well. 

The GA evaluates (14. ip in each evolutionary step for 
every chromosome of the population. The one with the 
best fitness will consequently have the smallest XsNia 
and will be the best candidate solution in its generation. 
Of all the steps in the execution of the algorithm, the 
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FIG. 2: \Bestfitness\ (lower is better) after 1000 generations 
as a function of the chromosome count in the population. 



evaluation of the fitness is the most expensive. GDF is 
appropriately modified to use (|4.ip as the basis of its fit- 
ness calculation. Additionally, GDF requires a number 
of parameters which determine the execution of the GA. 
These include the selection (percentage of the popula- 
tion that goes unchanged into the next generation) and 
mutation rates, the genome length of chromosomes in 
the population, and the total chromosome count in one 
generation. A maximum number of generations is also 
specified. We keep the default values for the selection 
rate (10%) and mutation rate (5%). 

The parameters which can significantly change the 
speed of the algorithm is the genome length and the 
population size. Given the use of GE, one should bear 
in mind that very long chromosomes with many genes 
are not always fully used in the transcription process to 
produce a formal expression. As such, we deem a chro- 
mosome length of 100 genes long enough, leading to solu- 
tions with a high goodness of fit, without being compu- 
tationally very expensive. Having fixed the chromosome 
length, we must choose an appropriate population size. 
Since the fitness function has to be evaluated for every 
chromosome in one generation, this is the most impor- 
tant parameter for the speed of execution. The GA is a 
stochastic algorithm and thus the speed of convergence 
to an optimal solution cannot be easily determined and 
is affected by a number of parameters. 

To find the most suitable population size, we make 
a preliminary scan of this parameter space for different 
number of chromosomes and determine what is the best 
fitness after 1000 generations. We determine sweet spots 
for 200, 410, 420, 440 and 480 chromosomes. The first 
of them is used as the chromosome count for the rest of 
the treatment, as it gives the fastest execution with good 
results. All tests use an evolution of 1000 generations. 

After the execution of the GA, we obtain an expression 
fioAiz) for the reduced distance modulus as the solution 
of best fitness and the corresponding XsNia- Using this 



7 




ACDM, I!,„=0.26 
GAfil 




0.0 



0.5 



1.0 



1.5 



FIG. 3: Left: The SNIa distance modulus against the redshift. The solid line represents the ACDM curve with Qrn ~ 0.26. 
The dashed line is a GA best fit function of eq. (|4.2p . Right: The data and the GA best-fit of eq. (|4.2I) residuals relative to 
ACDM. 



we can, through differentiation, obtain parameters of in- 
terest, such as dark energy equation of state parameters 
in the context of some modeL As already mentioned, 
this method gives no direct way to estimate the errors 
for the derived parameters. One cannot expect to get an 
estimate by just running the algorithm many times and 
obtaining slightly different parameters. The GA usually 
tends to converge at the same solution for a given dataset, 
unless we change significantly the population size or the 
number of generations. A way to circumvent this prob- 
lem is to use Monte Carlo simulation to produce synthetic 
datasets and rerun the algorithm on them. We can thus 
obtain a statistical sample of parameter values which will 
allow us to estimate the error. We sketch the procedure 
we follow: 

1. The GA is applied on the original SNIa dataset 
with the chosen execution parameters and a solu- 
tion for P,ga{z) is obtained. 

2. Choosing a particular dark energy model, we deter- 
mine its DE parameters from P-gaIz) using differ- 
entiations. These parameters will be used as seed 
values to produce synthetic datasets. 

3. A large number of synthetic datasets is produced 
from the DE model, using the seed values from the 
previous step. We add gaussian noise to the data 
points to simulate real observational data. The 
standard deviation for noise is taken to be equal 
to that of the corresponding observational point. 

4. The GA is rerun for the synthetic datasets. 

5. A new set of parameter values is determined from 
the solutions of the synthetic datasets. 

6. Outliers are eliminated using recursive trimming to 
reduce variance of the statistical sample. 

7. The final mean values and standard deviations of 
the DE parameters are determined. 



Using the above steps, we can obtain error estimates for 
the DE model we wish to examine. Step 6 is essential to 
reduce the number of outliers in the sample. The GA is 
stochastic and non-linear, so adding gaussian noise to the 
input data doesn't lead to a gaussian distribution of the 
output. For some of the synthetic datasets the algorithm 
will fail to converge fast enough, yielding a suboptimal 
solution, hence the presence of outliers. By recursively 
trimming the sample, we can significantly improve the 
accuracy of the result. 

B. Applications 

We will now give a number of examples where the 
method is applied and the results interpreted in the con- 
text of some dark energy model. We first use the GA to 
treat the more popular ACDM model, which is also an 
easier example to demonstrate the method. In the next 
section we consider a model with a dark energy equation 
of state parameter which depends on the redshift and 
determine the appropriate DE parameters of the model. 

1. ACDM 

For ACDM, the properties of dark energy are fully de- 
termined, since w = —1 and the only free parameter of 
the model is the total matter density il„i- As described 
above, we first run the algorithm on the SNIa dataset to 
obtain a seed solution for the reduced distance modulus, 
and thus the distance modulus ^ga{z)- Applying the GA 
to the original dataset as it is was found to produce solu- 
tions with unphysical behavior for small redshift, where 
we would expect H{z) Hq. For this reason, a fiducial 
point was added in the dataset at small z with a very 
low standard deviation, so that every good candidate so- 
lution should take it into account. The total during 
the GA runs included this artificial point. However, in 
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our subsequent treatment and final results we reverted 
to the initial dataset (all values mentioned later do 
not include the fiducial point). Using this method, the 
distance modulus was found to be 



[iGA{z) = log (z 



2.17145 



H'-«' + z + e/))+Mo. (4.2) 



Note that eq. (|4.2p is not unique since each time the GA 
is ran it may give a different but equally good and model 
independent approximation to the reduced distance mod- 
ulus /i(z). This function is then fitted against the corre- 
sponding modulus predicted by KCDM 



with 



IJ-KCDM [z) = 5 login Dl [z; n„i) + /xq , 



H {zf = (1 + zf + l- Qr^ 



(4.3) 



(4.4) 



to obtain the seed value of fim- We then use the Monte 
Carlo method to produce 100 synthetic datasets from 
(|4.3p adding gaussian noise and rerun the GA on them. 
The resulting set of fj,GA{z) from the synthetic datasets 
is used to obtain in the same fashion a statistical sample 
of rim values. Due to the non-linearity of the GA and GE 
implementation, the sample is not normally distributed 
and contains a lot of outliers. This is easy to establish 
both graphically through histograms, or by comparing 
the mean and standard deviation of the sample with its 
corresponding robust estimators, the median and median 
deviation. The large discrepancy between these quanti- 
ties is a sign of outliers in the dataset. We eliminate these 
by recursively removing outlier points. The heuristic we 
use is to discard in each step data points that are more 
than 2(7 away from the mean of the distribution at that 
stage. The trimming continues until we reach a station- 
ary state, where no other points are removed. From the 
improved statistical sample we determine the final values 
of the mean and standard deviation for the dark matter 
density to be flm — 0.308 ± 0.029, in good agreement 
with the WMAP value of f7,„ = 0.258 ± 0.030 Mi. 



TABLE I; Comparison of the mean and standard deviation 
of the Q.m statistical sample produced from the synthetic 
datasets and the corresponding median and median deviation. 
The initial large discrepancy is greatly reduced by trimming 
outlier points. 





Mean Median A (%) 


a 


Cmedian A (%) 


Initial 
Trimmed 


0.633 0.357 77.5 
0.308 0.298 3.26 


0.584 
0.029 


0.074 692 
0.016 83.3 



0.0 0.5 1.0 1.5 2.0 



00 0.2 0.4 0.6 0.8 1.0 



FIG. 4: Histograms showing the dispersion of the set of flm 
values. Upper: Initial set. Lower: After the recursive trim- 
ming of outliers. 



2. Chevallier-Polarski-Linder DE ansatz 



In this dark energy model, the equation of state pa- 
rameter is given by the ansatz. [43l. |4^ 



w (z) = Wo + wi 



1 + z 



(4.5) 



with two parameters wq and wi determining the temporal 
profile of w;(z). This model has been extensively used in 
previous parametric reconstructions of dark energy from 
SNIa data. The distance modulus of the model is given 
by equations (|2.ip and (|2.2p using 



H{zf = H^[nm{l + zf 

+ (1 - a„) (1 + z)3(i+-o+"'i) 



Here we use the matter density as an input parameter 
and attempt to reconstruct the dark energy profile by 
determining the two DE parameters wq and wi. We re- 
peat the process three times using as ftm the WMAP 
value ±1(7. As we see in Table 2, the WMAP value gives 
the best goodness of fit. The method is roughly the same 
as for KCDM , only now we have to determine two pa- 
rameters. After running the GA on the SNIa dataset 
and getting the seed values for wq and wi, we produce 
100 synthetic datasets and reapply the algorithm to get 
a sample of parameter values. Again, we have to remove 
outlier points to decrease the dispersion of the data. To 
do so, we recursively remove data points lying outside 
the 90% confidence level in every step, until we reach a 
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FIG. 5: The DE equation of state w{z) for different values of fim- 

TABLE II: Values for mean and median and their corresponding errors for the Chevallier-Polarski-Linder DE ansatz. Large 
discrepancies between the mean and the median signal the presence of outliers. Overall can be significantly -and adversely- 
affected by such points. The trimmed samples generally present a considerably better behavior. 







Mean wo 


Median wq 


Awo (%) Aawo (%) 


Mean wi 


Median wi 


Awi (%) A 


awi (%) 




Initial 


0.23 


-1.8 ±6.9 


-2.4 ± 1.3 


24.6 


424 


4 ±25 


5.5 ±5.5 


23.5 


345 


352.2 


Trimmed 


0.23 


-2.3 ±2.0 


-2.4 ± 1.1 


4.7 


73 


6.0 ±7.9 


5.6 ±3.8 


7.7 


109 


471.0 


Initial 


0.26 


-9 ±69 


2.4 ± 1.5 


278 


4480 


30 ± 240 


5.3 ±7.5 


456 


3140 


3797 


Trimmed 


0.26 


-1.8 ± 1.6 


-2.2 ± 1.1 


18.4 


37.8 


3.7 ±6.5 


5.2 ±4.7 


28.8 


37.8 


343.4 


Initial 


0.29 


-50 ± 470 


-2.4 ± 1.2 


1960 


38500 


200 ± 1700 


4.6 ± 6.6 


3810 


25500 


6667 


Trimmed 


0.29 


-2.3 ± 1.4 


-2.5 ±0.8 


7.9 


75.5 


4.4 ±5.5 


5.0 ±3.9 


12.1 


48.8 


420.6 



stationary state. From the improved sample values we 
get the mean and errors for wq and wi. 

The profiles we get for the DE equation of state pa- 
rameter are in accordance with previous parametric re- 
constructions based on this model. We see in Fig. 5 there 
is a trend for a redshift-varying wlz), with a crossing of 
the w = —1 line at late times. Still, the errors are large 
for the method to make any conclusive arguments about 
these effects. The errors are essentially due to the limited 
number of MC synthetic datasets used for this treatment. 
Increasing the number of datasets would lead to better 
statistics and decrease the errors further, but at a great 
computational cost, since the GA should then be applied 
to a much larger collection of data. 



of any parametrization in the treatment of the initial 
data which can reduce any model-dependent bias. In 
our case the desired function to be found is the reduced 
distance modulus It should be noted that running 

the GA several times may give different but equally good 
model independent approximations to the distance mod- 
ulus so what we actually need some statistics and 
a model to interpret the results. The reason for this is 
that if we want to constrain a parameter like Omegam, 
then already we have made a (minimal) assumption for 
the underlying theory. A similar approach has been also 
used in the literature before, see for example [32j| where 
the authors use the GA to discriminate between the var- 
ious SUSY models. 



V. CONCLUSIONS 

We presented a method for non-parametric reconstruc- 
tion of the dark energy equation of state parameter 
based on the genetic algorithms paradigm and gram- 
matical evolution. The method was applied to deter- 
mine parameters of cosmological models such as the 
AC DM and the Chevallier-Polarski-Linder DE ansatz. 
This approach produces results that are in agreement 
with previous parametric and non-parametric reconstruc- 
tions. The main advantage of the GA is the absence 



On the other hand, the shortcomings of the method in- 
clude computational speed issues, since the execution of 
the GA is usually very expensive regarding the CPU time 
needed. This in turn limits the amount of precision one 
can obtain using Monte Carlo methods to estimate errors 
for the final model parameters. However, this issue can 
soon be resolved with the advent of faster systems and 
more efficient GA implementations. Other open issues in- 
clude the determination of optimal algorithm parameters 
for the SNIa dataset, as well an extension of the method 
to analyze a larger variety of dark energy models. 
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